%% Initialization
%clear all
close all
clear global
clc
delete(gcp('nocreate'))
distcomp.feature( 'LocalUseMpiexec', false )

% restoredefaultpath;
rehash toolboxcache;

%% Add paths
addpath('../Matlab')

%% Read network (GRN+interface nodes)
G = importdata(strcat('adjmat-contextualized-',subfix,'-interface.txt'));
InData = importdata(strcat('pheno-contextualized-',subfix,'-interface.txt'));

%Formatting
RowName = G.textdata;
RowName(1) = [];
G = G.data;
BoolData = InData.data;
GNames = InData.textdata; % the same as RowName

initial_state = 1; % from experimental data

% Nodes indexes
delim = '\t';
to_perturb = importdata(strcat('to-perturb-',subfix,'.txt'),delim);
node_ind = zeros(size(to_perturb));
for i = 1:size(to_perturb)
    node_ind(i) = find(strcmp(GNames, to_perturb(i)));
end

% State of nodes to perturb
pert_BoolData = BoolData(node_ind, initial_state);
ups = node_ind(pert_BoolData==1); %index in BoolData(pert_ind, initial_state)
downs = node_ind(pert_BoolData==0);
both = node_ind(pert_BoolData==2);

% Unperturbed network (no nodes connected)
GRN_ind = ~ismember(1:size(G,1), both); %logical vector to select only GRN
unp_G = G(GRN_ind, GRN_ind);
unp_BoolData = BoolData(GRN_ind,:);

% Unperturbed attractor
[Att_states, count] = boolSim_pert(unp_G,unp_BoolData(:, initial_state),1000,[]);
m = [sum(Att_states==unp_BoolData(:, initial_state)) / size(Att_states, 1)*100, ... 
       sum(Att_states==unp_BoolData(:, initial_state)), count, sum(~(Att_states == unp_BoolData(:, initial_state)))];
unperturbed_Att_states = Att_states;
m_flipped = strjoin(GNames(find(Att_states ~= unp_BoolData(:, initial_state))),';');

%% Perturbations
% perturb the TFs already in the GRN
singles = ups;
Mups = zeros(size(singles,1),4);
Mups_flipped = cell(size(singles,1),1);

for i = 1:size(singles,1)
    ind_pert = singles(i,:);

    % Simulation
    [Att_states, count] = boolSim_pert(unp_G,unp_BoolData(:, initial_state),1000,ind_pert); 
    
    %Compare with unperturbed attractor
    expected_Att_states = unperturbed_Att_states;
    expected_Att_states(ind_pert) = ~expected_Att_states(ind_pert);
    flipped_genes = sum(~(Att_states == expected_Att_states));
    
    Mups(i,:) = [sum(Att_states==expected_Att_states) / size(Att_states, 1)*100, ... 
        sum(Att_states==expected_Att_states), count, flipped_genes]; 

    Mups_flipped{i} = strjoin(GNames(find(Att_states ~= expected_Att_states)),';');
end

singles = downs;
Mdowns = zeros(size(singles,1),4);
Mdowns_flipped = cell(size(singles,1),1);

for i = 1:size(singles,1)
    ind_pert = singles(i,:);

    % Simulation
    [Att_states, count] = boolSim_pert(unp_G,unp_BoolData(:, initial_state),1000,ind_pert); 
    
    %Compare with unperturbed attractor
    expected_Att_states = unperturbed_Att_states;
    expected_Att_states(ind_pert) = ~expected_Att_states(ind_pert);
    flipped_genes = sum(~(Att_states == expected_Att_states));
    
    Mdowns(i,:) = [sum(Att_states==expected_Att_states) / size(Att_states, 1)*100, ... 
        sum(Att_states==expected_Att_states), count, flipped_genes]; 

    Mdowns_flipped{i} = strjoin(GNames(find(Att_states ~= expected_Att_states)),';');
end

% add and perturb the interface nodes without known state
singles = both;
% set to 1
Mboth_up = zeros(size(singles,1),4);
Mboth_up_flipped = cell(size(singles,1),1);
for i = 1:size(singles,1)
    ind_pert = singles(i,:);
    pert_G = G(~ismember(1:size(G,1), both(~ismember(both,ind_pert))), ~ismember(1:size(G,1), both(~ismember(both,ind_pert))));
    pert_BoolData = BoolData(~ismember(1:size(G,1), both(~ismember(both,ind_pert))),:);
    pert_BoolData(pert_BoolData(:,initial_state) ==2,initial_state) = 1;
    
   % Simulation
    [Att_states, count] =boolSim_pert(pert_G,pert_BoolData(:, initial_state),1000,[]);
    
    %Compare with unperturbed attractor
    comparable_Att_states = Att_states(GRN_ind,:);
    flipped_genes = sum(~(comparable_Att_states == unperturbed_Att_states));
    
    Mboth_up(i,:) = [sum(comparable_Att_states==unperturbed_Att_states) / size(Att_states, 1)*100, ... 
        sum(comparable_Att_states==unperturbed_Att_states), count, flipped_genes]; 
    Mboth_up_flipped{i} = strjoin(GNames(find(comparable_Att_states ~= unperturbed_Att_states)),';');
end

% set to 0
Mboth_down = zeros(size(singles,1),4);
Mboth_down_flipped = cell(size(singles,1),1);
for i = 1:size(singles,1)
    ind_pert = singles(i,:);
    pert_G = G(~ismember(1:size(G,1), both(~ismember(both,ind_pert))), ~ismember(1:size(G,1), both(~ismember(both,ind_pert))));
    pert_BoolData = BoolData(~ismember(1:size(G,1), both(~ismember(both,ind_pert))),:);
    pert_BoolData(pert_BoolData(:,initial_state) ==2,initial_state) = 0;
    
    % Simulation
    [Att_states, count] =boolSim_pert(pert_G,pert_BoolData(:, initial_state),1000,[]);
    
    %Compare with unperturbed attractor
    comparable_Att_states = Att_states(GRN_ind,:);
    flipped_genes = sum(~(comparable_Att_states == unperturbed_Att_states));
    
    Mboth_down(i,:) = [sum(comparable_Att_states==unperturbed_Att_states) / size(Att_states, 1)*100, ... 
        sum(comparable_Att_states==unperturbed_Att_states), count, flipped_genes]; 
    Mboth_down_flipped{i} = strjoin(GNames(find(comparable_Att_states ~= unperturbed_Att_states)),';');
end

M = [Mups;Mdowns;Mboth_up;Mboth_down];
M_flipped = [Mups_flipped;Mdowns_flipped;Mboth_up_flipped;Mboth_down_flipped];
names = cat(1,strcat(GNames(ups),'_0'),strcat(GNames(downs),'_1'),strcat(GNames(both),'_1'),strcat(GNames(both),'_0'));

%% Export the results
fName = strcat(subfix,'-perturbation_1_interfaceNode.txt');  
fid = fopen(fName,'w');            
if fid ~= -1
    fprintf(fid,'%s\t%s\t%s\t%s\t%s\t%s\n','Perturbation', 'Match_fraction_initial', ...
        'Match_number_initial','Simulation_count','Flipped_genes','Names');
    for j=1:size(M,1)
        fprintf(fid,'%s\t',names{j});       
        fprintf(fid,'%3.2f\t%d\t%d\t%d\t',M(j,1), M(j,2), M(j,3), M(j,4));      
	fprintf(fid,'%s\n',M_flipped{j}); % print the list of flipped genes
    end
    % unperturbed
    fprintf(fid,'%s\t','Unperturbed');
    fprintf(fid,'%3.2f\t%d\t%d\t%d\t',m(1,1), m(1,2), m(1,3),m(1,4));     
    fprintf(fid,'%s\n',m_flipped); 
    %fprintf(fid,'\n');
    fclose(fid);     
end


